#include<iostream>
#include<vector>
#include<queue>
#include<malloc.h>
#include <math.h>

void V_cycle_one(double *f,double *u, int n, int ng,int ncycle_M ,double RA,int RO,int IO,double *au,double *y)
//f是右侧函数对应的向量,u是初始残差,n是区间个数加一,n=(2^ng)+1,ncycle_M是最大迭代次数,RA是相对精确度,RO是限制算子的类型,IO是插值算子的类型,au是程序计算出的最后的结果.


{
void addint_one_l(double *uf, double *uc, double *res, int nf);
void addint_one_q(double *uf, double *uc, double *res, int nf);
void copy_one(double *aout, double *ain, int n);
void fill0_one(double *u, int n);
void interp_one_l(double *uf, double *uc, int nf);
void interp_one_q(double *uf, double *uc, int nf);
void relax_one(double *u, double *rhs, int n);
void resid_one(double *res, double *u, double *rhs, int n);
 void rstrct_one_fw(double *uc, double *uf, int nc);
void rstrct_one_i(double *uc, double *uf, int nc);
void slvsml_one(double *u, double *rhs);
 double maxnorm_one(double a,double *u,int n);

 unsigned int jcycle,jj,jpost,jpre,nn,i;
 double Maxnorm,umaxnorm,a=-1,h,h2i,pumaxnorm,rateR;
 double *ires[16],*irhs[16],*iu[16],*mu,*er;

 //er为误差=真实值减去计算值

  h=1.0/(n-1);
 h2i=1.0/(h*h);
mu=(double *)malloc(sizeof(double)*2000);
 
er=(double *)malloc(sizeof(double)*2000);
 
 for (jcycle=1;jcycle<=ncycle_M;jcycle++) {
  std::cout<<"The "<<jcycle<<"-th V cycle:"<<std::endl;
  std::cout<<""<<std::endl;
  
  //ires,irhs,iu全部给空间
  for (jj=ng;jj>=1;jj--) {
    nn=pow(2.0,jj)+1;
    irhs[jj]=(double *)malloc(sizeof(double)*2000);
    ires[jj]=(double *)malloc(sizeof(double)*2000);//残差
    iu[jj]=(double *)malloc(sizeof(double)*2000);
    fill0_one(iu[jj],nn);fill0_one(ires[jj],nn);fill0_one(irhs[jj],nn);
  }
  
  for (jj=ng;jj>=2;jj--) {
    nn=pow(2.0,jj)+1;
     
    if (jj==ng){
      copy_one(irhs[jj],u,n);
      }
    
    for (jpre=1;jpre<=2;jpre++){
      relax_one(iu[jj],irhs[jj],nn);
      }
    
    resid_one(ires[jj],iu[jj],irhs[jj],nn);
    nn=(nn+1)/2;
  
    if (RO==0)
      rstrct_one_fw(irhs[jj-1],ires[jj],nn);
    else
      rstrct_one_i(irhs[jj-1],ires[jj],nn);
  }
  
  slvsml_one(iu[1],irhs[1]);
  nn=3;
  
  for (jj=2;jj<=ng;jj++) {
    nn=2*nn-1;
    if (IO==0)
      addint_one_l(iu[jj],iu[jj-1],ires[jj],nn);
    else
      addint_one_q(iu[jj],iu[jj-1],ires[jj],nn);
  
    for (jpost=1;jpost<=2;jpost++)
      relax_one(iu[jj],irhs[jj],nn);
     }
  
  //u是残差,pumaxnorm是残差的最大范数
   if (jcycle==1)
   pumaxnorm=maxnorm_one(a,u,n);
  else
    pumaxnorm=umaxnorm;
  //重新计算精确值和残差
  
  for (i=2;i<n;i++){
    au[i]=au[i]+iu[ng][i];
  }
  for (i=2;i<n;i++){
   u[i]=-h2i*(au[i-1]+au[i+1]-2*au[i])+f[i];
  }
   u[1]=-h2i*(au[2]-2*au[1])+f[1];
   u[n]=-h2i*(au[n-1]-2*au[n])+f[n];

   //输出残差的最大范数
   for (i=2;i<n;i++){
     mu[i]=u[i];
   }
   umaxnorm=maxnorm_one(a,mu,n);
   rateR=(pumaxnorm-umaxnorm)/pumaxnorm;
   std::cout<<"The maximum norm of residuals after "<<jcycle<<" times V cycle="<<umaxnorm<<std::endl;
   std::cout<<""<<std::endl;
   std::cout<<"The "<<jcycle<<"-th residual reduction rate is about "<<rateR<<std::endl;
   std::cout<<""<<std::endl;
  //判断精确度是否达到
  
   for (i=1;i<=n;i++){
     er[i]=y[i]-au[i];
   }
  Maxnorm=maxnorm_one(a,er,n);
  std::cout<<"The maximum norm of error after "<<jcycle<<" times V cycle="<<Maxnorm<<std::endl;
  std::cout<<""<<std::endl;
  if (umaxnorm<RA) break;
   
 }

}

